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Abstract 

In  this  note  we  present  an  implementation  of  the  fast  marching  algorithm  for  solving 
Eikonal  equations  that  reduces  the  original  run-time  from  0(1V  log  N)  to  linear.  This 
lower  run-time  cost  is  obtained  while  keeping  an  error  bound  of  the  same  order  of 
magnitude  as  the  original  algorithm.  This  improvement  is  achieved  introducing  the 
straight  forward  untidy  priority  queue ,  obtained  via  a  quantization  of  the  priorities 
in  the  marching  computation.  We  present  the  underlying  framework,  estimations 
on  the  error,  and  examples  showing  the  usefulness  of  the  proposed  approach. 

Key  words:  Fast  marching,  Hamilton  Jacobi  and  Eikonal  equations,  distance  func¬ 
tions,  untidy  priority  queue. 


1  Introduction 


The  fast  marching  method  [13]  has  been  introduced  to  solve  the  static  Hamilton- 
Jacobi  (Eikonal)  equation  \VT\F  =  1,  in  a  computationally  efficient  way.  Here 
F  >  0  is  the  front  moving  speed  and  T  is  the  travel  time.  F  and  T  may  also 
be  interpreted  as  travel  cost  and  intrinsic  distance  respectively.  Regarding  the 
computationally  efficient  implementation  of  this  equation,  Tsitsiklis  [16]  first 
described  an  optimal-control  approach,  while  independently  Sethian  [12]  and 
Helmscn  [5]  both  developed  techniques  based  on  upwind  numerical  schemes. 
The  complexity  of  their  approach  is  0(N  log  N),  where  N  is  the  total  number 
of  grid  points. 1  The  algorithm  is  an  extension  to  the  classical  Dijkstra  tech¬ 
nique,  and  is  based  on  finding  (at  each  iteration)  the  point  with  the  minimal 
T  value  in  the  narrow  band  set  of  points  that  are  being  updated,  and  setting 
it  to  be  alive  (points  that  already  got  their  definitive  T  value).  Neighbors  of 
this  point  are  then  updated  and  the  process  is  repeated  until  all  points  in  the 

1  In  practice,  the  number  of  visited  points  during  the  computation. 
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domain  have  been  processed.  All  involved  numerical  computations  use  upwind 
schemes.  The  point  with  the  minimal  T  value  is  usually  found  using  a  heap 
priority  queue  based  on  a  tree.  Insertion  time  to  such  a  heap  is  O ( log  n)  where 
n  is  the  number  of  points  in  the  narrow  band.  Due  to  the  positiveness  of  F, 
every  point  is  visited  at  most  a  constant  number  of  times,  and  from  this  the 
time  complexity  of  the  algorithm  is  0(N  log  N).  See  [13]  for  additional  details. 

Due  to  the  broad  applications  of  weighted  distance  functions  obtained  by 
solving  the  above  mentioned  Eikonal  equation,  e.g.,  [6,13],  its  efficient  im¬ 
plementation  and  numerous  extensions  have  been  studied.  In  particular,  fast 
sweeping  algorithms  have  been  proposed,  e.g.,  [2,15,17]  (early  ideas  in  this 
direction  were  proposed  by  Danielsson).  This  technique  is  based  on  using  a 
pre-defined  sweep  strategy,  replacing  the  use  of  the  heap  priority  queue  to 
find  the  next  point  to  process,  and  thereby  reducing  the  overall  complexity  to 
O(N).  2 

Here  we  propose  a  new  implementation  of  the  fast  marching  algorithm  which 
reduces  the  computational  complexity  to  O(N)  (being  in  practice  N  as  in  the 
original  fast  marching,  the  number  of  visited  points).  This  is  based  on  the  con¬ 
cept  that  when  solving  the  Eikonal  equation  for  the  current  grid  point,  we  can 
quantize  (round)  the  priority  values,  thereby  allowing  to  use  a  table  instead 
of  a  tree,  reducing  the  updating  complexity  from  O(logA)  to  0(1).  This  is 
done  with  the  help  of  a  data  structure  denoted  as  untidy  priority  queue.  We 
show  that  the  possible  error  introduced  by  this  simplification  can  be  kept  of 
the  same  order  of  magnitude  as  the  numerical  error  introduced  by  the  spatial 
discretization  inherent  to  numerical  implementations,  while  in  practice,  the 
errors  introduced  by  this  approximation  are  virtually  insignificant.  In  layman 
words,  the  idea  here  proposed  is  that  in  the  same  way  that  numerical  imple¬ 
mentations  introduce  errors  due  to  space  discretization,  we  can  also  allow  for 
errors  in  the  value  computed  from  the  Eikonal  equation  when  used  to  set  the 
priority  for  the  current  grid  point.  When  this  is  properly  done,  computational 
complexity  is  improved  at  no  theoretical  or  practical  cost.  The  rest  of  this 
note  provides  details  on  this  and  examples. 


2  Algorithm  Description 


Our  method  is  based  on  the  special  properties  that  grid  points  hold  in  the 
maintained  priority  queue.  When  the  T  value  of  all  points  in  the  queue  are 

2  Here  N  is  the  number  of  grid  points  in  a  bounding  box  of  the  region  of  interest. 
Although  theoretically  this  is  the  same  “A”  as  in  the  original  fast  marching  imple¬ 
mentation,  it  can  be  much  larger  in  practice,  in  particular  for  largely  anisotropic 
speeds  F  (computations  are  done  only  at  visited  points  in  fast  marching). 
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Fig.  1.  Operational  queue  for  the  narrow  band,  see  text  for  details. 

always  larger  or  equal  than  the  T  value  of  the  latest  point  extracted  from  the 
queue,  the  queue  is  denominated  as  a  monotone  priority  queue.  This  important 
subclass  of  priority  queues  is  used  in  applications  such  as  event  scheduling  and 
discrete  event  simulation.  In  the  fast  marching  algorithm  the  queue  maintained 
for  the  narrow  band  ( NB )  has  a  monotone  behavior.  If  we  assume  that  F  is 
bounded,  we  can  derive  that  the  T  values  of  the  points  in  the  queue  are  also 
bounded  since  there  is  a  maximal  possible  increment  Imax  over  the  point  with 
the  minimal  T  value  in  queue.  3  When  a  priority  queue  involves  only  points 
with  T  values  in  a  fixed  size  range,  it  is  possible  to  use  data  structures  based 
on  a  circular  array  (see  Figure  1).  Each  entry  of  the  circular  array  contains  a 
list  of  points  with  “similar”  T  value.  The  calender  queue  [3]  and  some  of  it’s 
improvements  [1,9,11]  are  all  such  queues  where  the  probabilistic  complexity 
of  the  insert  and  remove  operations  are  0(1).  4 

Let  Pmin  be  the  point  with  time  Tmin  which  is  the  minimal  time  in  the  queue 
and  Pr  with  time  Tr  be  the  next  point  to  be  removed  from  the  queue.  For  the 
calendar  queue  Pr  =  Pmin  always  holds,  i.e.,  we  are  always  guaranteed  to  get 
the  point  with  the  lowest  priority.  However,  since  the  fast  marching  algorithm 
already  has  an  error  dependent  on  the  grid  density  h,  0(h) ,  it  is  not  necessary 
to  strictly  keep  Pr  =  Pmin  to  achieve  the  same  general  order  of  numerical  error. 
We  propose  to  use  an  untidy  priority  queue  where  Pr  ~  Pmin-  The  queue  is 
based  on  the  circular  array,  which  is  a  simplification  of  the  calender  queue. 
Each  entry  in  the  circular  array  of  size  d  represents  a  discrete  (quantized)  level 
of  T,  uniformly  distributed.  Let  T  be  the  quantized  value  of  T  and  A  be  the 
difference  between  any  two  consecutive  discrete  levels.  Let  T0,  T\, ... ,  T<i- 1  be 
the  discrete  levels  of  the  array  entries  such  that  Tq  <  Ty  <  ...  <  Td-y.  The 
entry  Ty  keeps  a  FIFO  list  of  grid  points  with  T  values  in  the  range  [T),  Tl+l] 
that  quantize  to  T). 

Figure  1(a)  gives  an  example  of  the  proposed  untidy  priority  queue.  In  this 

3  The  value  of  Imax  depends  on  the  neighborhood  scheme  used,  e.g.  for  4-neighbors, 
Imax  —  max  (F). 

4  Yet,  worst  case  complexity  is  O(logn)  which  happens  if  too  many  elements  have 
nearly  equal  priorities. 
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example  d  =  6  so  that  6 A  =  Imax ■  Priority  values  T  are  then  assumed  to  be  in 
the  range  [T0,  T0  +  Imax\-  The  queue  keeps  track  of  the  memory  location  L0  of 
entry  T0  in  the  array  (the  entry  corresponding  to  the  current  value  of  T ).  In 
this  example,  Tq  is  the  3rd  entry  from  the  top  so  To  3.  When  a  new  point 
Pa  with  time  Ta  is  added  to  the  queue,  according  to  it’s  quantized  time  value 
Ta,  it  is  placed  in  the  end  of  the  list  of  entry  Tr  Where  j  ( Ta  —  T0)/A. 

Finding  entry  Tj  in  memory  can  be  done  using  a  simple  modulo  operation: 
Lj  =  ( j  +  L0)/A  mod  d.  The  next  point  to  be  removed  from  the  queue,  Pr, 
will  always  be  the  first  point  in  the  To  entry.  The  figure  displays  numbers 
which  correspond  to  the  order  in  which  the  points  would  be  removed  if  no 
new  points  are  added  to  the  queue  during  the  process.  Note  that  due  to  the 
quantization,  Pr  is  not  necessarily  Pmm,  but  we  can  guarantee  that  Tmin  =  Tr 
so  0  <  Tr  —  Tmm  <  A.  When  the  list  at  position  T0  is  emptied,  the  next  non¬ 
empty  list  is  used  (list  at  position  Ti  in  this  example).  The  entry  T0  is  now 
used  to  store  T  values  quantized  to  T0  +  6 A  (circular  queue).  For  consistency, 
all  labels  are  shifted  forward  one  position  so  the  new  Pmin  is  now  at  position 
T0  as  shown  in  Figure  1(b). 

Note  that  the  quantization  is  used  only  to  place  the  grid  point  in  the  queue, 
while  the  actual  T  value  is  used  to  solve  the  Eikonal  equation  when  the  grid 
point  is  selected.  Therefore,  errors  can  only  occur  due  to  wrong  selection  order 
(misorder),  see  below  for  more  details. 

The  average  complexity  of  the  remove  operation  is  0(1)  as  long  as  0(d)  < 
0(n)  (since  the  operation  may  involve  searching  for  a  non-empty  queue).  Se¬ 
lecting  a  constant  size  d  of  order  0(1)  or  by  using  automatic  resizing  techniques 
as  presented  in  [3],  it  is  possible  to  guarantee  a  worst  case  average  complex¬ 
ity  of  0(1).  The  insert  operation  has  no  searching  involved  and  therefore  it’s 
run-time  complexity  is  0(1). 

An  extreme  case  occurs  when  Imax  ^  0  and  increments  with  same  magnitude 
of  Imax  are  rare.  This  causes  a  “waste”  in  accuracy,  since  many  discrete  levels 
in  the  circular  array  would  contain  empty  lists.  Prior  knowledge  of  such  in¬ 
crement  distributions  permits  to  handle  the  rare  large  increments  separately 
as  suggested  in  [3] ,  thereby  avoiding  such  a  situation. 

To  conclude,  let  us  point  out  that  Tsitsiklis  also  described  0(N )  variations 
for  solving  the  Eikonal  equation  (although  the  constant  can  be  much  larger 
than  in  his  O(NlogN)  approach).  His  approach  is  based  on  buckets  and  a 
more  elaborated  discretization  of  the  equation,  which  is  “more  cumbersome 
and  is  unlikely  to  be  used  when  the  dimension  is  higher  than  three”  [16]. 
His  bucket  data  structure  is  a  linear  array  of  length  0(maxT/AT).  In  con¬ 
trast,  our  approach  does  not  require  a  new  discretization  (which  in  the  case 
of  [16]  includes  an  optimization  step  which  is  more  expensive  to  solve),  uses  a 
cyclic  array  thereby  no  needing  to  estimate  the  maximal  distance  and  being 
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more  memory  efficient,  has  very  small  constant  in  the  O(N)  complexity,  and 
quantizes  the  values  only  for  queue  placing  and  not  for  actual  computation. 
Thereby,  the  approach  here  proposed  is  simpler,  faster,  valid  for  any  dimension 
(as  well  as  for  the  extensions  reported  in  the  literature  for  computing  geodesics 
on  manifolds  [7,10]),  and  memory  efficient.  The  authors  of  [8]  also  mentioned 
that  numerical  precision  can  be  exploited  to  eliminate  the  tree  search,  but  no 
algorithm  was  proposed. 


2. 1  Error  Bounds 


Assuming  a  given  state  of  the  algorithm,  we  would  like  to  bound  the  additional 
error  in  the  computation  of  T  of  an  Alive  point  when  it  is  selected  out  of  order 
using  the  untidy  priority  queue. 

Let  Tro  be  the  minimal  possible  time  value  the  point  Pr  would  get  if  Pmin 
was  selected  (the  selection  an  accurate  queue/tree  would  give).  The  additional 
error  due  to  a  single  misordering  is  £  =  Tr—Tro.  Since  the  point  Pr  was  selected, 
the  untidy  priority  queue  restricts  Tr  and  Tnun  to  share  the  same  discrete  level. 
Hence,  Tr  —  Tmin  <  A.  Since  F  >  0,  it  is  possible  to  show  from  the  properties 
of  the  algorithm  that  Tmin  <  Tro.  When  combining  this,  we  obtain  that  the 
possible  error  clue  to  a  single  misordering  is  £  =  Tr  —  Tro  <  Tr  —  Tmm  <  A. 
Note  that  £  >  0  since  the  algorithm  may  only  decrease  the  value  of  points  in 
NB,  so  removing  a  point  early  may  only  cause  Tr  not  to  reach  Tro.  Therefore 
the  introduced  error  can  not  be  negative. 

Since  £  <  A  =  Imax/d,  by  selecting  d  =  0(1/ h)  we  can  achieve  e  =  0(h). 
This  is  the  same  order  of  magnitude  of  error  as  in  the  original  fast  marching 
algorithm,  only  that  computational  complexity  is  reduced  to  0(N ). 


3  Numerical  Experiments 


Numerical  experiments  show  that  errors  caused  due  to  the  misordering  are 
very  rare,  small,  and  in  practice  far  from  the  upper  bounds  computed  above. 
Figure  2  presents  error  statistics  versus  the  number  of  discrete  levels.  Numeri¬ 
cal  experiments  show  that  few  discrete  levels  are  enough  to  achieve  a  negligible 
error. 

Figure  3  gives  an  example  of  a  fast  marching  application.  The  task  is  to 
segment  the  lake  just  by  giving  4  points  on  the  shore  of  the  lake,  following  [4] . 
The  segmentation  based  on  the  distance  function  generated  using  the  untidy 
priority  queue  is  identical  to  the  segmentation  archived  using  an  accurate 
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(a)  (b) 

Fig.  2.  Error  statistics  compared  to  an  accurate  queue,  (a)  Statistics  when  creating 
a  time/distance  function  on  a  2000x2000  grid  with  random  0  <  F  <  1.  (b)  Same 
statistics  when  F  is  a  decreasing  function  of  the  image  gradient,  Figure  3. 


Fig.  3.  The  segmentation  example  shows  that  the  quantization  error  does  not  inter¬ 
fere  with  fundamental  fast  marching  applications.  (Image  size  1163  x  1463.) 

queue.  Yet,  the  run-time  when  using  the  untidy  priority  queue  is  much  shorter 
(see  below).  Additional  examples  are  available  online  at 
http:  /  /  mountains.ece.umn.edu/~liron/fastmarching/. 

Figure  4  shows  the  0(N )  run-time  complexity.  The  average  number  of  queue 
operations  done  per  point  in  the  gird  is  constant  as  the  grid  size  N  goes  to 
infinity.  Using  our  implementation,  the  time  required  to  compute  distances 
(from  a  single  seed)  in  a  2000x2000  grid  with  random  0  <  F  <  1  using  1000 
discrete  bins  is  1.251  seconds.  When  F  is  an  edge  map  obtained  from  the 
1163x1463  image  displayed  in  Figure  3(a),  the  time  required  is  0.591  seconds. 
The  computer  used  for  the  measurements  is  a  Pentium  4  laptop  with  2Ghz 
processor  speed. 

Acknowledgments-  We  thank  Gustavo  Brown  for  ideas  on  queue  based  irn- 
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Fig.  4.  Run-time  statistics.  Only  queue  search  operations  are  considered.  For  the 
heap  operations,  these  include  both  the  heap  access  and  the  priority  comparisons. 
For  the  untidy  priority  queue,  operations  are  array  entry  access.  Statistics  done  on 
a  square  grid  with  a  random  F  function. 

plementations  and  Ron  Kinnncl  for  reminding  us  of  the  O(N)  algorithm  in  [16]. 
This  work  was  partially  supported  by  the  Office  of  Naval  Research,  the  Na¬ 
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